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Abstract: We study the holographic field theory dual of a probe SU{2) Yang-Mills field in a 
background (4 + l)-dimensional asymptotically Anti-de Sitter space. We find a new ground state 
when a magnetic component of the gauge field is larger than a critical value. The ground state forms 
a triangular Abrikosov lattice in the spatial directions perpendicular to the magnetic field. The 
lattice is composed of superconducting vortices induced by the condensation of a charged vector 
operator. We perform this calculation both at finite temperature and at zero temperature with a 
hard wall cutoff dual to a confining gauge theory. The study of this state may be of relevance to 
both holographic condensed matter models as well as to heavy ion physics. The results shown here 
provide support for the proposal that such a ground state may be found in the QCD vacuum when 
a large magnetic field is present. 
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1 Introduction 

The study of black hole instabilities is an important research topic that has led to very interesting 
results. In particular, within gauge/gravity duality, the study of Anti-de Sitter black hole solutions 
and their stability properties is important for understanding thermal states on the gauge theory 
side. Last year in [1], some of the authors of the present paper studied an SU(2) Einstein- Yang- 
Mills model at finite temperature in asymptotically AdS space. They found that when a magnetic 
component of the gauge field reaches a critical value in units of the temperature, the system becomes 
unstable. Though the critical value of the magnetic field for onset of the instability was calculated, 
the new ground state of the system was not known. In the current work we calculate a ground 
state solution, using a perturbative analysis similar to the one performed by Abrikosov in [2] for 
type II superconductors. In agreement with the work of Abrikosov we find that the ground state 
is a triangular lattice. There have been many attempts recently to model lattices holographically 
with the goal of providing more realistic models for condensed matter systems [3, 4], and this 
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novel procedure for generating a lattice dynamically adds to these developments. Moreover, our 
holographic model provides support for recent QCD studies of p meson condensation from a strong 
magnetic field [5-7]. The effect described here is similar to the Nielsen-Olsen solution for gluon 
condensation [8] and to magnetically catalysed W boson condensation [9-11]. 

Although this is the first holographic calculation to explicitly uncover an Abrikosov lattice in 
3+1 dimensions, it is not the first to examine spatially inhomogeneous phases of strongly coupled 
field theories. In [3], the authors studied the holographic construction of an Einstein-Maxwell-scalar 
theory at finite temperature and density. They looked at the gauge theory optical conductivity, 
which is the conductivity in the direction of an applied electric field. They broke the translational 
invariance explicitly by imposing scalar field boundary conditions in the form of a lattice modulated 
in one of the Minkowski spatial directions. A fully backreacted solution was found which thus 
induces a spatially inhomogeneous black hole solution. This leads to an extremely rich behaviour 
of the frequency dependent optical conductivity. At low frequencies there appears a Drude peak. A 
Drude peak is a broadening of the zero frequency delta peak in the conductivity. In real materials 
this is due to impurities and finite temperature effects. The Drude peak is not present when 
translational invariance is unbroken. The solution also exhibits a power law behaviour at frequencies 
intermediate with respect to the temperature, and a constant value in the high frequency regime. 
The power law behaviour is the same as that found experimentally in cuprates, while the constant 
value at high frequencies is expected from conformal invariance. The setup in this context was a 
(2+l)-dimensional model where the lattice was periodic in only one of the spatial dimensions. A 
more realistic lattice structure would be highly desirable. 

While the lattice in the approach of [3] was implemented in the boundary conditions, there are 
a number of other mechanisms known that lead dynamically to ground states without translational 
symmetry. One approach was pioneered in [12] by studying a Yang-Mills-Chern-Simons theory in 
the gauge/gravity context. It was shown that the Chern-Simons term can induce an instability 
which leads to a ground state with both translational and rotational symmetry breaking. Such 
work was continued in [13-17] where the spatially homogeneous phase was found to be unstable 
in a variety of gravitational contexts in the presence of Chern-Simons couplings. The perturbative 
analysis of quasinormal modes that become tachyonic at finite momentum gives a relatively simple 
computational tool for finding instabilities to ground states without translational symmetry. These 
solutions are found to induce a helical current [18-20]. Interestingly, a Chern-Simons term is not 
always enough to induce such an instability. It was shown by [21] that this type of instability does 
not exist in the D3/D7 system. 

Translational invariance can also be broken with a magnetic field or with magnetic monopoles. 
The former was first studied by Gauntlett et al. in [22] and [23] where the instability of the 
magnetically charged black hole in a top-down framework was studied in detail. In the latter work, 
an infinite family of solutions coming from D = 11 supergravity was shown to exhibit a magnetically 
catalysed instability. Such work is important as it proves that these instabilities can also come from 
real string theory constructions. The subject of magnetic monopoles in (3 + l)-dimensional AdS 
space was studied in [24] and [25]. These magnetic monopoles are solutions to the scalar field in a 
Yang-Mills-Higgs theory with an SU (2) gauge field. In a certain limit where the monopole magnetic 
charge becomes large and a "monopole wall" is formed, it was shown in [24] that there is a W boson 
instability. In [25] a hexagonal lattice ground state of these monopole walls was found numerically. 

In [26] the (3 + l)-dimensional gravitational dual of a type II superconductor was constructed. 
A type II superconductor is one for which the external applied magnetic field has two critical 
values. When the magnitude of the magnetic field increases beyond the lower of the two critical 
values, the field starts to penetrate the superconducting condensate. Some of the condensate 
remains until the magnitude of the field is increased beyond the upper critical value, at which 
point superconductivity is completely destroyed. Just before the upper critical value is reached 
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from below, the ground state of the system is a triangular Abrikosov lattice [2]. The authors 
of [26] constructed a holographic superconductor with the same properties and found the Abrikosov 
lattice ground state explicitly. This is different from our model, however, because here we find a 
superconducting Abrikosov lattice ground state that is induced by an SU{2) magnetic field, rather 
than being destroyed by it. Moreover, in contrast to this model, we do not need a finite density. 
Our model is a cousin of holographic p-wave superconductors where the condensation is induced by 
a finite isospin density, holographically realised by a non-trivial temporal component of the SU(2) 
gauge field (see [27] and [28, 29] as well as the recent [30]). Here, in contrast, a spatial component 
of this gauge field has a non-trivial profile. Whereas in [28, 29], a Meissner effect in shown to occur 
by which a magnetic field reduces the transition temperature, here it is again the magnetic field 
which induces condensation at zero density. 

In addition to being interesting in the broader context of holographic lattices, the model we 
discuss serves as supporting evidence for a phenomenon first described by Chernodub et al. in [5, 6]. 
There it was proposed that the QCD x QED vacuum may itself be susceptible to a superconducting 
transition when a magnetic field of the order of the QCD scale is present. Such extreme conditions 
are rare but they may be present for a few femtoseconds during highly off-centre heavy ion collisions. 
The discovery of this phase came about through the study of an effective field theory description 
(the DSGS model proposed by Djukanovic, Schindlcr, Gegelia and Scherer in [31]]) of p mesons 
interacting with a magnetic field. A destabilisation of the vacuum was shown that would clearly 
lead to the condensation of charged and neutral p mesons. This breaks the U(l) gauge symmetry 
and leads to a superconductor with the quark-antiquark pairs in the mesons acting as Cooper pairs. 
The instability was also found using an extended Nambu-Jona-Lasinio model with SU(3) colour 
and SU(2) flavour in [32]. Lattice gauge theory studies were then performed looking at QCD in 
strong magnetic fields and these indicate the same instability. Moreover, using the DSGS model and 
guided by the Ginzburg-Landau model of type II superconductors, a solution was found in which the 
p meson condensate forms an Abrikosov lattice made up of superconducting vortices [7] . This may 
be relevant experimentally. Evidence has mounted at both RHIC and the ALICE experiment at 
CERN that strong magnetic fields may contribute to the physics of the strongly coupled quark gluon 
plasma as charges are quickly accelerated during the interaction period [33, 34]. The importance 
of these effects remains a contested topic because the time scales involved are small. However, 
given that strong magnetic fields may be present, it is interesting to ask if traces of this p meson 
condensate could be detected. 

In the current paper we find a possible ground state of the system in [I]. As mentioned above, 
this system was shown to be unstable under the imposition of a large SU(2) magnetic field. 1 We 
show that it has very similar properties to the ground state of a type II superconductor near the 
upper critical magnetic field as well as to the ground state in the model of Chernodub et al. In 
other words, the ground state is a triangular Abrikosov lattice. We take here a very simple model 
of a strongly coupled finite temperature quantum field theory in (3 + I )-dimensions with a global 
SU(2) symmetry. The dual gravity theory is an SU(2) Einstein- Yang-Mills theory in (4 + I)- 
dimensions with a magnetic component of the SU (2) switched on. We work entirely within the 
probe approximation, which means that the Yang-Mills term is small compared to the Einstein- 
Hilbert term in the action. We also fix the gauge in such a way that the gauge theory condensate is 
transformed under a U(l) subgroup of the global SU(2) symmetry. There appears to be a certain 
universality to the triangular lattice ground state. Here we show that it forms in both the AdS 
Schwarzschild background (dual to a finite temperature field theory) as well as the hard wall cutoff 
model (dual to a confining field theory). It would be interesting to uncover exactly how universal 

It was shown in [35] that the same sort of instability occurs in the Sakai-Sugimoto model, but there the ground 
state has also not been found. 
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these results are. 

The two holographic models that we study have several important differences from QCD. In 
the hnitc temperature model there is no confinement or chiral symmetry breaking and so there are 
no goldstone bosons (pions) present which are the normal decay modes of the p meson in QCD. The 
hard wall model has its conformal symmetry broken only by an IR boundary condition which sets 
a confinement scale. However, the phenomenology of these two models appears to be close enough 
to that of QCD to compare qualitatively with the models of Chernodub et al. 

In section 2 we provide the details of the holographic setup. There we also outline the per- 
turbative analysis of fluctuations of the SU(2) gauge field near the critical magnetic field. Since 
we follow the philosophy of Abrikosov's calculation of the ground state in type II superconductors, 
which was done in the Ginzburg-Landau model, in section 3 we give a brief outline of this approach 
and then follow it to solve perturbatively up to third order. In section 4 we discuss the numerical 
results and analyse the free energy of the different lattice solutions, showing that the triangular 
lattice has the lowest free energy of all the Abrikosov solutions studied. It is important to note that 
we are not able to show conclusively that we have found the ground state but we are able to find 
a state with lower free energy than the translationally invariant state and that has lowest energy 
within a large class of lattice solutions. In section 5 we conclude and give an outline of important 
future work. 

2 Holographic setup 

2.1 The finite temperature and hard wall backgrounds 

The system we study is an Einstein- Yang-Mills theory on the (Poincare patch of) an asymptotically 
AdSs geometry with an SU(2) gauge field. The action is 



where g is the Yang-Mills coupling, Gat is the 5D gravitational constant and L is the AdSs radius. 
R and F are the Ricci scalar and Yang-Mills field strength respectively. 

We consider the probe approximation, where the Yang-Mills term is small compared to the 
Einstcin-Hilbert term, so that the gauge fields do not backreact on the geometry. We thus choose 
a fixed 5-dimensional background metric, given by 



where the asymptotically AdS region is at u — > 0. We study two different models. The first is a 
finite temperature model where the background is AdS Schwarzschild, first proposed in [36]. In 

4 

this case, f(u) = 1 — -V, where uh is the location of the planar black hole horizon. The Hawking 
temperature of the black hole is T — 1/ttuh- The second model is the hard wall cutoff model, 
proposed in [37, 38], where f(u) = 1 and the geometry terminates at a radial distance uc- This 
model corresponds to a zero temperature theory (ujj = oo), but it still has a scale tip which 
corresponds to a confinement scale in the gauge theory. The intrinsic scales in these theories allow 
us to form a dimensionless magnitude for the magnetic field. This will be the parameter that we 
tune in order to find the instability of the spatially invariant ground state. Without loss of generality 
we can choose units where uh = 1 in the finite temperature theory and uc = 1 in the confining 
theory. Factors of uh and uq can then be restored through dimensional analysis. In the following 
the exact form of the metric is not important until we come to solving the numerical equations in 
the radial direction of AdS. 




(2.1) 




(2.2) 
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S = -j^ I d 5 x^9 tr (F^Fn , (2.3) 



2.2 The Yang-Mills action 

The relevant part of the action simplifies to 

1 

with the equations of motion 

VF£ V + e abc A b,l F° v = . (2.4) 
The SU(2) gauge field is A = AtT a dx^, for a = 1 . . .3. We use the convention where the Lie 



algebra basis is given by r" = with a a the Pauli matrices, and the structure constants f abc 
are defined by [r a ,r fc ] = ie abc T c so that f abc = e abc . With these definitions, the components of the 
field-strength tensor F = dA + A A A become 

F% = V*2 - duAl + e ahc A\Al . (2.5) 

It will be important to understand how gauge transformations affect the system. Under a gauge 
transformation e ( x A transforms as 

A^A^ + SA^ = e iK A»e- ik - id^e^ . (2.6) 

When Kix 11 ) is an infinitesimal transformation, this becomes 

SA% = V^A a = d^A a + e abc AlA c . (2.7) 

The gauge transformations give us the freedom to fix the gauge A^ = 0. We work in this gauge 
from now on. 

In this paper we look at the effect of a strong (flavour-)magnetic field given by F^ y = B, with 
all other components of F^ vanishing. As we will see, when B becomes large 2 , other components 
of F become non-zero dynamically. To get a consistent set of equations we therefore consider a 
gauge field A of the form 

A= A°(x,y,u)T a dx» . (2.8) 

It turns out that we can turn off the t and z dependence of the gauge field and still have consistent 
equations. This simplifies the equations. Turning off the t dependence guarantees a static solution. 
Turning off the z dependence, where the z direction is parallel to the magnetic field, yields a lattice 
in the x, y-plane. 

2.3 Perturbative expansion of the gauge fields 

The SU(2) gauge freedom in the choice of A needs to be handled carefully, as discussed in [1]. 
There the authors take infinitesimal fluctuations of A and form linear combinations that transform 
covariantly in the fundamental of a U(l) subgroup of SU{2), and then fix the remaining gauge 
symmetry. The [/(l)-covariant fields condense and thus, on breaking the U(l) symmetry, this leads 
to a superconducting 3 ground state. The [/(l)-covariant fields, in the gauge where A = xBr 3 dy 
and F® = for the non-superconducting solution, are 

E L = x 2 BE x + i{xd x E y -E y ) , (2.9) 



2 Since we have chosen the units where ujj = 1 or uq = 1, B is a dimensionless quantity. Restoring the units, the 
statement is that Bv? H = B/(ttT) 2 or Bv? c ~ ^I^-qcd ^ s l ar g e i or that B is large compared to the radial scale of 
the background. 

3 Note that this is actually a superfluid because the gauge symmetry is mapped to a global symmetry on the 
boundary. 
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and its complex conjugate El- Here = A 1 + iA 2 ^ are small fluctuations of the gauge fields. We 
will however not use the El in the calculations because it is technically easier to work with E^, 
but the resulting solution can be used to calculated E L - See [1] for more details. 

Substituting equation 2.8 into equation 2.4 yields nine coupled partial differential equations in 
the variables x, y and u. Of these nine equations of motion, six are dynamical equations for each 
field Ax^' 3 , and three equations are constraints. The constraint equations arise from the equations 
of motion for the components A 1 ^ 2 ' 3 , which were chosen to be zero using gauge symmetry. 

In solving the PDE's, we follow the strategy of [2, 39], which works as follows. When the 
magnetic field B is smaller than some critical value B c , the field configuration A 3 = xB, A 3 = 
and A].' 2 = solves the equations of motion. This is the normal phase of the superconductor. 
As shown in [1], the system enters a new phase when the magnetic field is increased beyond some 
critical value B c . In this phase, the ground state has a non-trivial profile for all fields in the 
ansatz equation 2.8. It is thus a superconducting phase because the non-trivial profile is dual to 
a condensate in the field theory that breaks the U(l) symmetry. We look for this configuration at 
some value of B infinitesimally above B c , where the condensate is still small and the equations can 
be expanded perturbatively in a parameter proportional to the overall condensate size. For more 
about this technique, see [1] and references therein. 

We are interested in the phenomenology of this system close to the phase transition, where B is 
infinitesimally above B c . In this case the correction to the trivial background solution will be small 
and a perturbative analysis can be made in a parameter e ~ (B — B c ). For notational convenience 
we leave this parameter s explicit when studying the expansion. However, it will be absorbed into 
the definition of the perturbative corrections to the fields when we come to minimising the energy. 
We thus write an ansatz for the expansion in the form 

A 3 y = xB c + eA 3 y +e 2 a 3 y + ..., (2.10) 
A^=eA^+e 2 a^ + ... for (a, M ) + (3, y) , (2.11) 

and solve the equations order by order in e, as detailed in section 3. 
2.4 Gauge field boundary conditions 

According to the holographic dictionary the gauge fields in the UV can be expanded in powers of 
u. From this expansion the expectation values and couplings of their field theory duals can be read 
off. The field A 3 y has a boundary expansion given by 

<»Uo = ^+<}« 2 + -" , (2-12) 

(3) 

where s x , y is the value of the source, which in this case is the externally applied magnetic field 

' (3) 

potential. t4,# is proportional to the expectation value corresponding to the magnetisation. We 
will set the boundary conditions so that the applied magnetic field is not corrected by the higher 
order perturbations in e, whereas the magnetisation will obtain a non-zero value. 
Similarly, the fields A\' 2 y have a boundary expansion given by 

I u -$.q = s x,y v x,y u +■•■ ' (2-13) 

fl 2) 

where s x ,' y corresponds to the source of the operator that will condense to break the U (1) symmetry. 
We adjust the boundary conditions in such a way that this vanishes so that the symmetry breaking 

(12).. 

is spontaneous. v x , y is proportional to the vacuum expectation value of this operator, which we 
will read off to find the resulting supercurrent in the superconducting phase. 

Boundary conditions are also imposed on the fields in the IR. In the case of the black hole 
background, we impose regularity at the horizon and in the case of the hard wall model we impose 
Neumann boundary conditions. 
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2.5 The gauge theory ground state energy 

In finding the ground state, it is important to be able to calculate the energy of the field theory 
solution from the action. We will want to compare the solutions in the normal phase to those 
in the superconducting phase. The energy T of the gauge theory solution is found by using the 
holographic dictionary. In the case of the finite temperature solution, we are in the canonical 
ensemble and we calculate the free energy, which is T = — Tin 2 = — S c i with our conventions. 
Here S c i = — ^-p / d b Xyf—gF^ u F a ^ 1 ' is the classical action. In the hard wall case, we are simply 
calculating the energy of the field configuration, which is defined in terms of the classical action in 
the same way. Since we are only interested in whether the energy of a particular superconducting 
solution is lower than that of the normal phase solution, we can simply calculate the difference 
AJ 7 = F s — J- n and thus do not need to implement holographic renormalisation. Here F s is the 
energy of the superconducting phase, while T n is the normal phase energy with A y = xB and all 
other components zero. We also need to take care of the fact that S c i diverges when we perform 
the integral over the Minkowski directions. This is easy to fix by considering the energy density 
f2, which is obtained by integrating S c i only over the world volume of one lattice cell and dividing 
by its volume. Having explained how to calculate the energy of a field configuration, in the next 
section we turn to the problem of solving the equations of motion to find the ground state. 

3 Solving the equations 

3.1 The comparison with Ginzburg-Landau theory 

Before turning to the equations of motion, it helps to first look at the Ginzburg-Landau equations 
for an analogy. In some suitable units defined in [2, 39], they are 

(-iv - iy ^ - *p + \^\ 2 ip = o , (3.i) 

V x V x A = -i^V^i-ipV^) - \?p\ 2 A . (3.2) 

Only the structure of these equations is important, so we have ignored constant factors. Here ip 
is the wave function of Cooper pairs, and A is the electromagnetic vector potential. The nine 
equations of motion in our system can be split into two groups that roughly correspond to the two 
equations above. 

The first of the two groups, hereafter called the condensate equations, contains the six equations 
for the fields A x ' 2 y u . The superconducting condensate of the dual field theory, which is like %j) above, 
is found by looking at the subleading terms in the boundary expansion of A x 'y, equation 2.13. Of 
the six equations in this group, the dynamical equations are for A x ,2 y and the constraint equations 
are for A]/ 2 . So this first group is analogous to equation 3.1. The analogy can be made more 
clear. As mentioned above, we should be solving for the field of equation 2.9, but it is technically 
easier to work with the field definitions £ x , y = A xy + iA xy . Making this definition, we are able 
to combine the six real equations into three complex equations, two dynamical and one constraint. 
The constraint equation relates £ x and £ y such that there is only one complex degree of freedom 
left, which is analogous to the state ip. All this is hard to see at the non-perturbative level, but 
it illustrates the strategy we follow for solving the equations at each order: we use the constraint 
equation to reduce the two dynamical equations into one, and then solve it. 

The second group of equations, which we call the magnetic field equations, is for the fields 
A x y u , corresponding to A in equation 3.2 above. There are three such equations, one of which is 
a constraint. At each order we will be able to use the constraint to separate the equations into one 
for A x and one for A y . 
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3.2 The gauge field perturbative expansion in more detail 

Having denned the ansatz for our gauge potential in equation 2.10 we can learn more about the per- 
turbative expansion by studying the non-linear structure of the equations of motion. The equation 
for Af, is 







(3.3) 



We see that the magnetic field components appear in the linear terms, while the condensate compo- 
nents appear in quadratic terms. This suggests that a contribution to the condensate components 
that is first order in the perturbative expansion influences a second order contribution in the mag- 
netic field components. More generally, an odd order contribution to the condensate components 
influences an even order contribution to the magnetic field components. 

This structure is common throughout all the equations of motion. It turns out that terms in 
the perturbative expansion of the magnetic field components that have an odd order vanish. The 
even order terms in the condensate components can then also be set to zero. We can thus constrain 
the expansion ansatz of equation 2.10 to 



eE x 



A 
A 



xB c + s 2 a y + 0(e /l ) , 



(3.4) 



Here the calligraphic letters denote the non-perturbative fields. E x . v and e x _ y are first and third 
order contributions to the condensate components, respectively, while a x y are second order correc- 
tions to A x y . 

Because of this convenient expansion of the fields, the condensate components and the magnetic 
components decouple at each order. That means that at each order, we only need to work with 
fields we have already solved at previous orders. Our strategy is thus to solve for the fields in the 
following sequence: 

£x,y = sE X y -\- £ 6i,j + O(e^) , 



Al = xB c 



A 



3 




■ 0{e A ) 
+ 0{e 4 ) 



In the next section we start with the linear order solution, which will shed more light on the 
procedure that must be implemented at higher orders. 



3.3 Solving the equations to linear order 

Using the expansion 3.4 and keeping terms to linear order, we find that there are six remaining 
equations given (in complex form) by 



= -iB c xd u E v - dyd u E y 



= B 2 x 2 E x - iB c E v 



- A d u E x - fd 2 u E x - 2iB c xd y E x 



dyE x + iB c xd x E y + d x d y E y , 



= 2iB c E x 



f' d U Ey - fdlEy + lB C X8 X E X + d X 8yE X - d^Ey 



(3.5) 

(3.6) 
(3.7) 
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Here, as above, f(u) — 1 — u A for the AdS Schwarzschild model and f(u) = 1 for the hard wall 
model. 

We can solve these equations by following Abrikosov [39]. The solution is given by 



E v = -iE, 



X 1 



(3.8) 



oc 



E x — 



]T C n e 



in, 



ky-lB e (x-^) 2 jj^ _ 



(3.9) 



n——oo 



U(u) is determined by solving 



u" + 



( 



/'(«) 



/(«) 




/(«) 



£7 = , 



(3.10) 



subject to the constraints 17(0) = and f7'(l) = 0. For the AdS Schwarzschild model (f(u) = 1— it 4 ), 
the latter constraint comes from imposing regularity at the horizon. It is possible to calculate B c 
by numerically finding the value at which U(u) satisfies these constraints. There is an infinite tower 
of solutions to B c , but we are only interested in the lowest one, which is where the phase transition 
occurs. For further details on solving this equation in the AdS Schwarzschild model, see [1]. For the 
hard wall model, /(it) = 1 so the equation simplifies to the extent that it can be solved analytically, 
the solution being a Bessel function. Qualitatively the solutions for U (u) in both models look very 
similar and we are only interested in their numerical form. For the AdS Schwarzschild model, we 
get B c ss 5.1, while we get B c sa 5.8 4 for the hard wall model. 

It should be noted that the solution 3.9 for E x agrees precisely (except for the factor of U(u)) 
with the linear order solution for the order parameter close to the upper critical magnetic field H C 2 
in the theory of type II superconductors, as seen in [39]. It is also the result found by Chcrnodub 
et al. in [7]. Depending on the values of the parameters C n and k (to be determined by the 
higher order equations in the perturbative expansion), E x corresponds to different inhomogeneous 
functions in the x, y-plane. We are particularly interested in finding those with lattice symmetries 
that represent evenly spaced vortices running in the z direction in the gauge theory. 

3.4 The Abrikosov lattice solution 

Before going beyond linear order, we discuss the possible solutions we can expect. The number of 
coefficients specifying a configuration can make the problem of finding the lowest energy solution 
unmanageable without making use of some symmetries. We can argue that, since nothing in the 
setup is explicitly breaking translational invariance in the x, y-directions, the solution should be a 
highly symmetric lattice. A nice review of how lattices can be formed from the Abrikosov solution 
(3.9) is given in [40]. There the authors explain that in order for \E X \ to be a lattice solution, the 
coefficients C n must have the same magnitude |C n | and moreover be periodic in some integer P, 
that is, C n — C n +p. 

In [2], Abrikosov first studied the simplest solution, a square lattice. In this case, P = 1, 
implying that C n = C for all n, and k — \/2ttB c . Later Kleiner et al. in [41] generalised the 
analysis by looking at P = 2, with C\ — ±iCo = ±zC. This choice of coefficients specify a general 
rhombic lattice, with the shape of the rhombus controlled by varying k. In particular, a square 
lattice can be obtained by choosing k — \/irB c . This square lattice is the same as Abrikosov's 
solution with P = 1, but it is rotated by tt/4 and translated. A triangular lattice is obtained by 
choosing k = 3^^/ttBc. 

4 This is the zero of the Bessel function of the first kind Jq(\/B). 
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Figure 1. A lattice cell, illustrating the meanings of L x and L y for a fixed area cell. 



To show how this works, we first substitute P = 2 and C\ = iCo = iC into the solution for E, 
which simplifies to 



It is then easy to see the symmetries \E x (x + [m + \q]L x ,y + [n + \q]L y )\ = \E x (x, y)\ for integers 
to, n and q. L x and L y are the lengths of the lattice cell in the x and y directions, and are given 
by L x = 2k/B c and L y = 2ir/k. See figure 1. 

We follow the approach of Kleiner et al, which is to compute the free energy density of the 
lattice for a range of values of the ratio L x /L y = k 2 /irB c . This essentially means that we vary 
k. The energy is computed numerically from the analytic expressions we obtain at each order 
in the following sections. What we find agrees with their result that the triangular lattice has 
the lowest free energy of the P < 2 solutions. When doing this, magnetic flux conservation is an 
important constraint. The total applied magnetic field per unit area is constant, and each lattice cell 
corresponds to a vortex with a single quantum of magnetic flux. This means that when comparing 
the energy of different lattices, we should make sure that they have the same magnetic flux per unit 
area, which in turn means that their lattice cells have the same area. Fortunately with this ansatz 
that is always the case since the area L x L y — 4n/B c is independent of k. 

In the following sections we calculate analytic expressions for the higher order corrections to 
the gauge field. We keep P and the coefficients C n general, except for imposing the periodicity 
condition C n = C„+p. 

3.5 Higher order contributions to the energy 

In order to find the ground state solution we must calculate the energy of the superconducting 
solutions and compare them to the normal phase solution. We can study the form of the energy as 
defined in section 2.5 to see how far we must go in the perturbative expansion of the gauge fields. 
The energy has terms that are quadratic and quartic in the gauge potential. The quartic term 
ensures that the energy is bounded below, because it has a positive coefficient. The quartic terms 
have lowest perturbative contributions of order e 4 . One might expect contributions of order e 3 



oc 




(3-11) 



TL— — OO 



coming from the zeroth order magnetic field contribution multiplied by three first order corrections. 
However, from equation 3.4 it can be shown that such terms do not arise. Thus we should expect to 
expand to third order in A x ,2 y and fourth order in A 3 y . However, it turns out that going to fourth 
order is not necessary because inserting ansatz 3.4 into the action of equation 2.3, we find that the 
only fourth order terms from A 3 y that appear at the fourth order of the action are proportional to 
~ dyd^ 3 — dxOy^ . Here a^ 3 and a y 4 ^ 3 are the fourth order corrections to A 3 , and A 3 , respectively. 
This term respects the lattice symmetries, thus on performing the integration over the lattice cell 
to get the free energy density, it vanishes by Stokes' theorem. 

We saw above that the parameters k and C n in the solution 3.9 are not fixed by the equations 
of motion to linear order. This is due to the fact that to linear order, the different vortices do not 
interact. We can therefore not expect to fix any of the coefficients C n or the spacing parameter k 
at this order. In fact, trying to calculate AQ to this order, which has no quartic terms in A, one 
finds that the free energy density is not bounded below; increasing the overall magnitude of the 
condensate always decreases Afi. To see which configuration, that is, which set of values for C„ 
and fc, is energetically favourable, we clearly have to go beyond linear order. 

3.6 Solving the equations to higher orders 

In this section we solve the equations of motion up to third order in the perturbation parameter. 

The second order corrections to the gauge fields contribute to the potentials A 3 and A 3 , that 
is, a x and a y in 3.4. These fields source the external magnetic field and the magnetisation. We 
impose that these corrections must vanish at the AdS boundary, so that the dual field theory has a 
constant applied magnetic field. We find however that they do not vanish throughout the bulk. In 
particular they develop non-vanishing subleading terms in the boundary expansion, representing a 
magnetisation in the field theory. 

In appendix A we explain how the equations for the Fourier modes of the fields a x and a y can 
be decoupled. This yields the following equations 

n ff a -3 / A f,2 2 4B c 2 mV S 

ud u [-d u a x y (m,n,u) ] - [ k n 



,22 




where 



and 



+ T^e- V e^Qd+n )U 2 = , (3.12) 



T x = -i ^^ n , T y = 2i n m , (3.13) 



a 3 x Jx, y, U )=£E e^ 1 "'* fi* iV (m, n, u) . (3.14) 



As before, P defines the periodicity in the C n . The parameters m and n correspond to the Fourier 
space levels of these fields. In order to calculate the solution a 3 y (x, y, u) we will in theory need to 
solve these equations for all values of m and n. However, it will turn out to be sufficient to only 
study the first few Fourier modes. The numerical procedure for solving these will be explained in 
section 3.7 

At third order we are studying the perturbative corrections to the condensate. Here we calculate 
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the corrections e x and e v . It is reasonable to assume that the answer is of the form 



oo 



sE x + e 3 e. 



X 



e J2 {C n U{u)+E?c x , n {u))e- inky -^*-^) 



(3.15) 



n— — oo 



oo 



sE v + e 3 e. 



v 



e J2 (-iC n U(u) + e 2 c y , n (u))e- ink ^ B °(*-^) 



(3.16) 



n— — oo 



where we have made use of equation 3.8 to relate the first order terms C n U(u) in £ x and £ y . 
c (x,y),n( u ) is the first perturbative correction to the condensate where the u dependence is a function 
of n in contrast to the first order term. 

We can write e x and e y in Fourier space, then use the three condensate equations discussed in 
section 3 to calculate these corrections. The one constraint equation can be used to decouple the 
other two equations. We then have one equation for c x ^ n {u) and one for c y ^ n {u). Further details 
are provided in appendix B. 

3.7 Numerical solutions 

Having separated the equations into ordinary differential equations in u by the method outlined in 
the appendices, we can now solve them numerically. Both the second and third order equations 
take the same general form, given by 



This equation can be solved numerically by picking some parameters for C n and k that give a 
particular lattice and then using a shooting method to integrate from u — 1 (the horizon/hard wall 
cutoff) to u = (the AdS boundary). It is an inhomogeneous second order differential equation, 
so there are two integration constants. The first is fixed by imposing regularity at the horizon or 
Neumann boundary conditions at the hard wall cutoff. This fixes the value of d u <p(l). The second 
constant is obtained by demanding that <f>(0) — 0, so that the fields vanish at the AdS boundary. 
This vanishing corresponds to both the magnetic field strength corrections and the source for the 
condensate being set to zero. We fulfil this boundary condition by adjusting </>(l). Unlike in the 
case of the first order equations, the equations here are not homogeneous and thus the source sets a 
scale with which the value </>(l) can be compared. Changing 4>(\) in this case thus acts as more than 
just a scaling for the solution and so is used as the tuning parameter to satisfy the UV constraint. 

For all of the equations, we can implement this procedure for arbitrary integers m and n, 
corresponding to the different Fourier modes of the gauge fields. This will then give a Fourier 
coefficient a x y (m, n, u) that can be used to determine a xy (x, y, u). Fortunately we do not have to 
do the calculation for many different values of m and n, because as the values get large, the source 
term gets suppressed exponentially. This can be seen in equation 3.12 for the second order terms 
and is true also for the third order equation. For a vanishing source, the equations for a x or d y 
have only the trivial solution. This means that a x y (m,n,u) is negligibly small for large m or n, 
and we can therefore truncate the Fourier series for a x y beyond m, n w 3. 

4 Results 

4.1 Finding the minimum energy state 

As explained above, we wish to find the values of the parameters k, P and C n = C n+ p that give 
the minimum energy state. These parameters define the shape of the lattice. Our analysis is only 
valid for B slightly above B c , where B c was determined in section 3.3. The first step is thus to 




(3.17) 
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pick a value for B in this vicinity. We then choose a set of lattice parameters that give us the 
lattice solution we wish to consider. As mentioned in [40], for lattice solutions all the C n must 
have the same magnitude C. We can therefore fix C n up to the normalisation C, along with a 
value of k, according to the discussion in section 3.4. We then substitute these values into the 
energy density that was defined in section 2.5. It takes the form Af2 = a\eC + a,2£ 2 C 2 + . . . . 
At this point we see that we can redefine C by absorbing a factor of e, which we call C £ . C £ is 
the only parameter left unfixed up to this point in the analysis. Here the 0.$ are values that are 
calculated numerically from substituting the solutions to the equations of motion into the expression 
for the energy derived in appendix C. AS1 forms a Mexican hat potential, which is easy to minimise 
numerically. An illustration of this procedure is shown in figure 2. The plot in figure 3 shows the 

AO. 
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Figure 2. The change in energy density as a function of C e , the overall condensate scale. The leftmost 
curve corresponds to B = B c , which is never negative for nonzero condensate. Curves for B < B c are 
similar. Increasing B beyond B c yields the curves to the right, and we see a clear minimum of the free 
energy that is lower than the free energy of the normal phase forming. The dashed line traces out the 
minimum of each of these curves, which corresponds to the energetically preferred size of the condensate 
as a function of B. This plot was generated in the AdS Schwarzschild model, but the hard wall model has 
qualitatively similar results. 

energy-minimising value of C £ as a function of magnetic field near the phase transition at B c . It 
shows that C e ~ [B — B c )i, so the condensate 5 has a critical exponent of 1/2. 

Having minimised with respect to C e for a given value of B and a given lattice configuration, 
we can plot the difference in the energy between the normal and superconducting states. Figure 4 
shows Afi, the difference between the free energy density in the superconducting and normal phases, 
as a function of external magnetic field for two different lattices. The first lattice is square, and the 
second is triangular. Both are described in section 3.4. 

The curves in figure 4 are the result of calculations in the AdS Schwarzschild model, but we get 
the same results up to a rescaling of the axes for the hard wall model. In the AdS Schwarzschild 
model, the critical magnetic field B c s» 5.1, while in the hard wall model B c ps 5.8. Each curve shows 
that the free energy density is proportional to (B — B c ) . This shows that the phase transition is 
second order, as expected if one looks at the analogous case in Ginzburg-Landau theory. There one 
can show ([42]) that the free energy is proportional to (T — T c ) 2 , where T c is the phase transition 
critical temperature. 

5 Note that only the combination sC is physically relevant, not C or e independently. 
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Figure 3. C £ ~ the overall condensate size for the AdS Schwarzschild solution in units of the temperature, 
as a function of the external magnetic field B. For B < B c , the condensate is zero, and for B slightly above 
B c , we see a (B — _B C )5 scaling behaviour. The plot for the hard wall model is the same, up to a scaling of 
the B and C E axes. 




Figure 4. The change in energy density (compared to the normal phase) for the triangular and square 
lattices as the external applied magnetic field is varied. The phase transition happens at B c w 5.1, which 
is where the coordinate axes are centred. Af2 squaro — Afi t riangic is so small that the two plots are almost 
on top of each other. This is for the AdS-Schwarzschild model, but the plots for the hard wall model are 
identical except for the scale on the axes. In the hard wall model, B c « 5.8. 

4.2 An analysis of P — 2 solutions 

We now specialise to the case where the periodicity of the C n is P = 2. This describes a general 
rhombic lattice solution which includes both the triangular and square lattices. The P = 1 square 
lattice can be found within the P = 2 solutions up to translation and rotation. We here perform 
the analysis done in [41] as described at the end of section 3.3. 

By looking at the form of equation 3.11, it is possible to see that the triangular lattice occurs 
for R = L x /L y — y3 and R = l/v3- In general, R and 1/R give the same lattice but with 
the x and y directions flipped. This is why figure 5 displays the symmetry ACl(R) = Aft(l/R). 
The triangular lattice corresponds to a global minimum of the energy as a function of i?, as seen 
from the figure. There is a local maximum for the square lattice, which is when R = 1. As 
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R — > oo (or R — > 0), the free energy increases. Intuitively one can understand this by making use 
of the properties of Abrikosov vortices that we understand from type II superconductors. These 
vortices repel. Since R — > oo and R — > correspond to elongating the rhombic lattice cell (while 
keeping the area constant) neighbouring vortices are squeezed together, and since they repel, this 
is energetically unfavourable. The energy difference as a function of R is plotted in figure 5. In 




Figure 5. The change in free energy density as a function of R = L x /L y , the ratio of side lengths of a 
constant area lattice cell. This plot is for the AdS Schwarzschild model, but the plot for the hard wall 
model is identical up to a rescaling of the axes. When R = 1, the lattice is square and the free energy 
achieves a local maximum. When R — \/3 and l/\/3, the lattice is triangular and the free energy is at a 
global minimum. Note that the plot has the symmetry Att(R) = AS1(1/_R), which simply corresponds to 
swapping the x, j/-axes. 

figure 6 we present the contour plot of the modulus squared of the condensate in the x, y-plane 
for the minimum energy solution corresponding to the triangular lattice. This is calculated by 
studying the boundary expansion of the field E x and reading off the normalisable term. We could 




Figure 6. A contour plot of the modulus squared of the field theory condensate dual to E x in the ground 
state triangular lattice. Darker colours mean smaller values. 

also plot the magnetisation of the ground state, which is found from the normalisable term in the 
boundary value expansion of d x a^ — d y a x . However, it takes the same form as the condensate and 
the numerics indicate that it differs only up to a scale. 
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5 Conclusion 



In this work we have found a likely ground state for the black hole Yang-Mills instability analysed 
in [1]. The solution, being of a lattice form, clearly has much potential for analysis in condensed 
matter models, where the breaking of translational invariance has already been shown to be very 
important in getting realistic phenomenology. As we have explained, it also has possible implications 
for heavy ion collider physics. There are a number of interesting areas where we could apply similar 
techniques and perform further calculations to elucidate the phenomenology of the ground state 
found here. 

It is certainly important to understand exactly how universal this result is. We have seen that 
the difference between the normal phase and superconducting phase energy density of this solution 
is, up to a scale, independent of the background geometry in the two models that we have studied 
here. It would be very useful to understand precisely where this universality stems from and find 
how much we can deform the gravity solutions until the Abrikosov lattice is no longer the ground 
state. 

In the present work we have analysed lattices with P = 1 and P = 2, corresponding to 
square and rhombic forms, respectively. Going to P = 3 requires a large increase in computational 
power. While this would be an interesting further calculation, the analogous cases of type II 
superconductivity and the model of Chernodub et al point to the triangular lattice being the true 
ground state. We thus expect higher P lattices to be energetically disfavoured. 

Having found this solution, there are some extensions that can be made. It will now be possible 
to study time dependent fluctuations about the ground state. In order to do this we would have to 
introduce a second perturbative parameter in addition to the parameter e used in the current work. 
This would be analogous to the parameter a' in a D-brane construction. This would allow us to 
study the transport properties of the lattice ground state, by looking at current-current correlation 
functions. If we wish to study the effect of the lattice on the shear viscosity to entropy density 
ratio, we would have to introduce gravitational back reaction in our model. Clearly this will be 
a much more involved calculation, with many non-linear couplings, but if we want to study the 
theory in a more realistic scenario, where the stress energy tensor also has a lattice structure, such 
a calculation would clearly be important. 

It is expected that if the QCD vacuum is unstable to p meson condensation in extremely off- 
centre heavy ion collisions, then the timescale of the instability would not be enough to form a 
well-defined lattice. Abrikosov vortices may form, but the magnetic field would likely drop below 
the critical value before they had time to arrange themselves into a lattice. It would be very 
interesting to perform a real-time calculation in order study the formation of the vortices and their 
movements as the magnetic field increased and decreased through the lifetime of a single off-centre 
collision. 
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A Deriving the equations for a? x 

We substitute the ansatz 3.4 into the full equations of motion and neglect terms beyond quadratic 
order in e. To get rid of all appearances of E y , we use the relation that E y = —iE x from 3.8. Then 
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we find that there are only three equations in which the fluctuations al y appear. We focus on those 
three equations. 

The simplest equation of the three is the constraint equation, which came from the equation of 
motion for A„. To quadratic order, this equation is simply 

d u d x al + d u d y a 3 y = . (A.l) 

The first thing to do is integrate by u. This gives an integration constant, but by the fact that both 
al and a 3 must vanish at u = 0, this integration constant vanishes. So the even simpler constraint 
is 

d x al + d y a 3 y = . (A.2) 
This is all we need to decouple the other two equations in al and a 3 . These equations now become 

= 3 -E x U%E x + 3 -E x U%E x - \iE x U 2 8 x E x + \iE x U*8 x E x 

+ ud u (J-d u al^ + dial + 814 , (A.3) 

= - B c xE x E x U 2 - l -iE x U 2 d y E x + ^iE x U 2 d y E x 

- *E x U 2 d x E x - *E x U 2 d x E x + ud u [^dua^j + 9^ + d 2 x a 3 y , (A.4) 

which are partial differential equations with sources that come from the linear order solutions. 
These two equations only differ by their source terms, so we will focus on a^. a y should be similar. 
Using the expression 3.9, we can see that the source term is periodic with y ~ y + 3? . a? x must 
have the same periodicity, so we can write it as a Fourier series, 



al 



[x,y,u)= J2 e- inky a 3 x (x,n,u) . (A.5) 



The equation becomes 

Y^^e-^(-^ +x ) 2 - k ^(- k -^ +x ) 2 knC m C n+m U 2 



k 2 n 2 a 3 x + ud u ( -d u a 3 x ) + d 2 x a 3 x = . (A.6) 



We notice that the source term in this equation is periodic in the a;-direction; x ~ x + This 
lets us expand Fourier series in x as well: 



al = ^2e i?s ^ x al(m,n,u) . (A. 7) 



Writing the source term as a series lets us then obtain the equation 3.12 for the coefficients 
al(m,n,u). 

Calling the source term S(x), the naive way of finding its Fourier coefficients is to use the 
formula 



B r — 



S n = Ykj Q e l ^*S(x) . (A.8) 

However, the source terms contains Gaussians, and those are much easier to integrate when the 
domain of integration is the entire real line. So we do the following trick. Doing a continuous 
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Fourier transform on a periodic function gives a sum of 5-functions, 

f dx e ipx S(x) = f dx e tpx ^ e~ i ^^ x S n 

The coefficients in front of the ^-functions are what we are looking for. We get 

I dx e lpx S(x) = -J^^ie~^ + ^ +1 ^7-^knC m C m+n U 2 . (A.fO) 
•* V c m n 

Using 

oo oo P— 1 

E /H= E E/( pm +o ( A - n ) 

■m— — oo m=— oo /— 

and then using the symmetry Ci+p — Ci, the only m- dependence remaining in the sum comes from 

ikPmp 

e B c . Making use of the identity 

oo oo 

J2 e lmq = 2Tt J2 5{q-2Tnn) (A.12) 



m— — oo m— — oo 



and 8(ax) — gives us the sum over <5-functions from A. 9. Then we can simply read off the 
coefficients S n . This gives us the equation 3.12. 

B Deriving the equations for c xn , c y>n 
The third order equations of motion are 

= ia 3 x E x U' + a 3 y E x U' - iE x Ud u a 3 x - E x Ud u a 3 y + iB c xd u e y + d y d u e y + d x d u e x , (B.l) 
= -iB c xa 3 x E x U - 2B c xa 3 y E x U - E X E 2 X U 3 - a 3 x U8 y E x + 2m 3 y Ud y E x - a 3 y Ud x E x 

- 2E x lId y a 3 x + iE x Ud y a 3 y + E x Ud x a 3 y + iB c e y - iB c xd x e y - d x d y e y 

- B 2 c x 2 e x + 2iB c xd y e x + d 2 e x + ud u ({^e^ (B.2) 

= B c xa 3 x E x U + iE x E 2 x U 3 - ia 3 x Ud y E x + 2a 3 x Ud x E x - ia 3 y Ud x E x 

+ iE x Ud y a 3 + E x Ud x a 3 - 2iE x lId x a 3 - 2iB c e x - iB c xd x e x - d x d y e x + d 2 x e y + ud u f~^« e 3/j ■ 

(B.3) 

The first of these is the constraint equation. We use it to relate e x and e y . In order to do this, we 
first simplify it by noticing that, since 



n= — oo 



we have that iB c xd u e y + d y d u e y = —id x d u e y . We can then integrate the entire equation by 
u, imposing vanishing boundary conditions at the AdS boundary. The constraint equation then 
simplifies to 

= -2iE x J x - 2E x J y + ia 3 x E x U + a 3 E x U + d x e x - id x e y , (B.5) 
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where 



J Xi y(x,y,u) = / U{u)dua xy (x,y,u)du 



(B.6) 



This allows us to eliminate e x in equation B.3 (after differentiating it by x). We write each function 
as a Fourier series in y and find an equation for the coefficients c y ^ n . At this point the equation 
still has an x dependence, which can be eliminated by multiplying the equation by (nk — B c x) to 
make it an even function in x and then integrating dx. In doing so we use the solution for E x 
and the form for e y given by B.4, as well as the Fourier series representation of the other functions. 
Once this is done, we are left with an equation for e y in the form 3.f7. 
The resulting equation for c y ^ n is 



E 



27r g (ifc 2 J=(n-r) + i3 c 7r<,) 



(-2 (k 2 Pr + 2iB c Trq) J x<qi 
kP 



C n - r ( (2ik 2 Pr - 4:B c nq) J y , qr + (2iB c Trqa 3 x + (-ik 2 Pr + 4B c tt 9 ) af) U 



kP 

(3B C + 2k 2 g{-2r + q)) C n+q C n+r C n - r+q U 3 



where 



B C C y n -\- 1id U I U Cy n 



U^daOi (q,r,u)du , 



for i = x,y. 

A similar procedure gives the constraint equation in terms of the coefficients, 



= Cx t n(u) - ic y , n {u) 



E 



PkB r 



ik 2 Pr + 2B c nq) C„_ 



(B.7) 



(B- 



x (2 J x , qr - 2iJ y>qr - (a 3 Xtqr ~ ia a ytqr )U(u)) j . (B.9) 
Once the coefficients c y _ n are found, we use this to calculate c x> „. 
C Calculating the energy 

The difference between the energy of the superconducting phase and that of the normal phase is 



AT= wJ d5x ^~ 9 ( F ^ F ° 



I superconducting 



normal 



(C.I) 



We calculate the energy density by averaging over the domain 0<y<^,0<x<j^,0<u<l 
and t,z € M. Since the integrand is independent of t and z, the averaging amounts to simply 
dropping the integration over those variables. In the following expression we use 



£x,y — A x v + i-A-x ,y ~~ ^ ] ^(x,y),n( u ) e " 2 B " ! 



(C.2) 
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we write A% = and Ay = xB 4- a^, and call the averaged energy AO. The result is 
4g 2 AQ = / du<Qi(u)-\- ^l^ijn^n^ u) + ^3 (m, n,ii) + 04(771, n,i£)] 

^ m,n— — 00 

00 ^ 

X! n 5 (m,ra,g,r,it) > , 



(C.3) 



m,n,p,q— — oo 



where 



°i =t^— yi-r 71 (fd u Cj,id u C jt i +CjyCjy) +3(^,^3;,/ -iC x ,iCy,i) 



kPu ^ 2 , 

i=0 \i=a:,i/ 



n 2 =- 



7-3/ \ 2i?m7T „ 

kna^(m, n, u) a, y \tri, n, u) 



y j \\d u a?(m,n, u)\ 



(C.4) 
(C.5) 



1 — P— 1 

fi a = \ P 4B 



2k 2 P 2 u 



E 

1=0 



i= 2 -p 2 ( (3k 2 mP + 2iBnir) a x (n, m, u)C X} i+ m C yi i 



+ a x (n, — m, u)C Vi i ((3k 2 mP + 2iBmr) C x j+ m — 2ik 2 mPC y ^+ m ) 
+ a y (n, —to, u)C x j+ m (—AiBnnCxj + (ik 2 mP + 6-Bmr) C y j) 

+ a y (n, m, u)C x j +m Cy } i (-ik 2 mP - SBnir)^ , 

1 [ttB fc2 ( m2 +" 2 ) 



(C.6) 



4fcPw V 2 
p-i 

^ ^ {Cy ,l-\-mCy J+7iC x ^iCxJ+m+n ^Cx^+mCy^+nCx^+m+nCy.l C x jC x J-]-m+nCy .l-\-mCy ,£-fn) 
(=0 



(C.7) 



/ — P-i 

„ V7T-B fc 2 m 2 i(2i + m)r»ir B„ 

f2 K = — 2^ e 4B p ~^ 



Pfcu 



/=n 
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In these expressions, C x . n and C y ^ n are functions of u. Their complex conjugates are given by C x . n 
and C y ^ n , respectively. All the infinite sums in the energy C.3 can be terminated at a small finite 
value because of exponential suppression in the fix 5 terms. 
In deriving this, it helps to make use of the formulae 
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where the latter is valid whenever h(x, m, n) = h(x + L,m + P,n + P) 
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